

setwd("/Users/italolopez/Dropbox/NEP Evaluation - The role of beliefs/JPE Final Version/JPE Submission/Replication package/Outputs/Tables")
############################################################################
############################################################################
### DV: zptevir_irtscore2 ###
############################################################################
############################################################################

################# MV: Home score ###################
set.seed(7688);
require(MASS)
a=0.085
b=0.243
acov <- matrix(c(
  0.002116, 0,
  0,0.000441
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = FALSE) 
M=mean(ab)
med=median(ab)

################# MV: Nurturing Index ###################
set.seed(7688);
require(MASS)
a=0.078
b=0.222
acov <- matrix(c(
  0.002116, 0,
  0,0.000441
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = FALSE) 
M=mean(ab)
med=median(ab)

################# MV: Positive Discipline ###################
set.seed(7688);
require(MASS)
a=0.100
b=0.159
acov <- matrix(c(
  0.002209, 0,
  0,0.000441
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = FALSE) 
M=mean(ab)
med=median(ab)

################# MV: Self-efficacy ###################
set.seed(7688);
require(MASS)
a=0.101
b=0.251
acov <- matrix(c(
  0.002304, 0,
  0,0.0004
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = FALSE) 
M=mean(ab)
med=median(ab)

################# MV: Pacotis ###################
set.seed(7688);
require(MASS)
a=0.104
b=0.130
acov <- matrix(c(
  0.002209, 0,
  0,0.000441
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = FALSE) 
M=mean(ab)
med=median(ab)

################# MV: Social Support ###################
set.seed(7688);
require(MASS)
a=0.083
b=0.108
acov <- matrix(c(
  0.002209, 0,
  0,0.0004
),2,2)
rep=20000
conf=95
pest=c(a,b)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = FALSE) 
M=mean(ab)
med=median(ab)

################# MV: Behavioral and Beliefs Index ###################
set.seed(7688);
a1=0.085 #estimated coefficient a1
a2=0.078 #estimated coefficient a2
a3=0.100 #estimated coefficient a1
a4=0.101 #estimated coefficient a2
a5=0.104 #estimated coefficient a1
a6=0.083 #estimated coefficient a2

b1=0.154 #estimated coefficient b1
b2=0.102 #estimated coefficient b2
b3=0.057 #estimated coefficient b2
b4=0.158 #estimated coefficient b2
b5=0.035 #estimated coefficient b2
b6=0.015 #estimated coefficient b2

a1std=0.046 #SE of coefficient a1
a2std=0.046 #SE of coefficient a2
a3std=0.047 #SE of coefficient a1
a4std=0.048 #SE of coefficient a2
a5std=0.047 #SE of coefficient a1
a6std=0.047 #SE of coefficient a2

b1std=0.022 #SE of coefficient b1
b2std=0.023 #SE of coefficient b2
b3std=0.021 #SE of coefficient b1
b4std=0.022 #SE of coefficient b2
b5std=0.021 #SE of coefficient b1
b6std=0.020 #SE of coefficient b2

rep=20000 #number of simulated values
conf=95 #confidence level
a1vec=rnorm(rep)*a1std+a1 #create vector of simulated a1 coefficients
a2vec=rnorm(rep)*a2std+a2 #create vector of simulated a2 coefficients
a3vec=rnorm(rep)*a3std+a3 #create vector of simulated a1 coefficients
a4vec=rnorm(rep)*a4std+a4 #create vector of simulated a2 coefficients
a5vec=rnorm(rep)*a5std+a5 #create vector of simulated a1 coefficients
a6vec=rnorm(rep)*a6std+a6 #create vector of simulated a2 coefficients

b1vec=rnorm(rep)*b1std+b1 #create vector of simulated b1 coefficients
b2vec=rnorm(rep)*b2std+b2 #create vector of simulated b2 coefficients
b3vec=rnorm(rep)*b3std+b3 #create vector of simulated b1 coefficients
b4vec=rnorm(rep)*b4std+b4 #create vector of simulated b2 coefficients
b5vec=rnorm(rep)*b5std+b5 #create vector of simulated b1 coefficients
b6vec=rnorm(rep)*b6std+b6 #create vector of simulated b2 coefficients

total=(a1vec*b1vec)+(a2vec*b2vec)+(a3vec*b3vec)+(a4vec*b4vec)+(a5vec*b5vec)+(a6vec*b6vec)# simulated total indirect effects
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(total,low) #lower limit of confidence interval
UL=quantile(total,upp) #upper quantile for confidence interval
LL_perc=(LL/0.094)*100
UL_perc=(UL/0.094)*100
LL4=format(LL_perc,digits=4)
UL4=format(UL_perc,digits=4)
print(paste("T10-CIs","LB%",LL4,"UB%",UL4 ,sep=";" ))
out <-c("T10-CIs",LL4,UL4)
capture.output(out, file = "T10-CIs.txt", append = TRUE) 
M=mean(ab)
med=median(ab)

